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From  2004  -  Present: 


•  USNA  has  begun  studies  specific  to  the  Chesapeake  Bay 

•  Current  efforts  span  five  departments: 

•  Physics,  Math,  Chemistry,  Oceanography  and  Naval  Architecture 

•  Efforts  have  begun  to  join  CBOS 

•  Instrumentation  has  begun  in  the  Severn  River  and  College  Creek 

•  Study  of  small  estuaries  -  “feeder”  systems  -  to  the  Bay  started 

•  Three  Trident  scholars:  Gillary  and  Brasher  (2004),  Boe  (2006) 

•  Differing  approaches  taken:  Normal  Mode  Analysis, 

•  Navier-Stokes  integration,  COMSOL  model 

•  100  modes  calculated  for  Dirichlet  and  Neumann  (2005) 

•  Feasibility  study  for  Galerkin  method  to  extract  £)(2006) 

•  Initial  Value  Problem/  Dual  Time  Problem  /  Dual  Position  Problem 

•  ~  10  monitoring  stations  needed  for  ~it(xv  t ) 

•  Concerns  about  proper  mesh  for  Bay 

•  How  to  handle  varying  wave  speed  in  the  Bay? 
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Fig.  1  The  Black  Sea  compliments  of  NASA  World  Wind. 


Fig.  2  Monterey  Bay  from  NASA  World  Wind.  The  open  boundary 
with  the  Pacific  clearly  visible. 


Fig.  3  Chesapeake  Bay  from  NASA  World  Wind.  The  Atlantic 
ocean  open  at  the  southern  end,  allowing  in  salt  water.  The  north 
end  dominated  by  fresh  water. 


Briefly,  the  formulation  leading  to  the  calculation  of 
fluid  flow  stems  from  the  realization  that  the  vector  fields 
can  be  derived  from  two  scalar  fields,  which  are  the  solu¬ 
tions  to  the  Helmhotz  equation  under  Dirichlet  and  Neu¬ 
mann  boundary  conditions  [4]. 

Tt  =  V  x  [(«¥')  + Vx(/i®)].  (1) 


and  5  illustrate  the  varying  levels  of  complexity  that  have 
been  solved  using  NMA. 

The  basic  unit  of  calculation  used  throughout  this  pa¬ 
per  is  the  normal  mode.  Like  the  modes  of  a  guitar  string 
or  an  organ  pipe,  systems  obeying  the  Helmholtz  equation 
and  Dirichlet  or  Neumann  boundary  conditions  will  res¬ 
onate  in  states  referred  as  ’’normal  modes”.  For  the  Chesa¬ 
peake  Bay,  the  modes  calculated  are  energy  potentials  whose 
gradients  and  curls  of  gradients  correspond  to  the  vector 
current  fields  found  in  fluid  mechanics  (it). 


Here  is  the  stream  potential  where, 

N  f-dY  d'F\ 

«r,  =  MD  =  (-w,-S;),  (2) 

and  ®  is  the  velocity  potential  where, 

_>  ,  ,  (d®  d®\ 

“"  =  (“-v)»  =  <3) 

with  (m,v)  representing  the  surface  current  velocities  in 
the  x  and  y  directions  respectively.  The  total  velocity  field 
is  composed: 


image  Processing  of 


the  Cfj 


Galerkin  Method: 


oo 


n— 0 


u(x,t)  =  E  a(t)nu(x)D  +  b(t)nu(x) 


'n  \  J  N,n 


a{t)m  =  fu(x,t)datu(x)Dmdn 


OO 

E 

n= 0  L 

OO 

E 

n= 0  L 

OO 


a(t)m  =  /  £  [a(t)nu(x)D  n  +  b(t)nu(x)Nn J 
a(i)m  =  £  a(t)nluDuDdQ,  +  b(t)nluNuDdQ, 


n=0 

a(^)m  =  ^nmfl(4' 


—  fu{x,  t)datau(x)DndQ 
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Partial  Galerkin  Method: 


d(t)m 

j^u(lt,t)u(lt)DndVL 

r  °c  f  i 

d(t)m 

n=0 

d(t)m 

oo  f  ,  r  -| 

H  a(t)„  dQ  +  6(t)„  kli.r  dfl 

n  y  Jn  JVt  D,n  D,m  K  Jn  JVt  N,n  D,m 

n= 0 

°C  r  n 

(i) 

d(t)m 

crn  a(£)_  +  /3n  &(£) 

D,nm  \  Jn  '  D,nm  V  Jn 

n= 0 

Note  that  the  a  n  and  tfn  exhibit  wavelet-like  responses  in  the  spectral 

D,nm  1  D,nm  1  1 

domain. 
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Dual  Time  Problem  (Initial  Value  Problem) 


F(x^  t)  o,n 

F(X)  t)  data 


f(x)n9(t)n 

oc  r  i  r 

E_  \An9(t)Dtn  +  Bn9(t)N>n]  [Cnf(X)D,n  +  Dnf(X)N,n\ 


OO 

E 

n= 0  L 

(X) 


n=0 


F(x,  t)data  E  ACng Dnf Dn  +  BCngNnf Dn  +  ADngDnf Nn  +  BDngNnf 


n*yN,nJ  D,n  n^D^rid  N,n 


n*y  N,nJ  N,n  _ 


—  f  F(x,tl)dataf{x)D^md9l 


OO 


n=0 


f  ^2  ACngDnfD,n  n9 N.nf D.n  DmUf N.n  n9N.nf [ 


n*J  N,nJ  D,n  1  ^  nX D,nJ  N,n 


nX N,nJ  N,n 


«(*i)n  =  ^Cnj(ti)Bin  +  BC„j(ii) 


'JV,n 


—  j  F{xAi)dataf{x)Dmd9l 
fd(ti)m  =  f  F(x,ti)dataf(x)Nrnd9l 
l(h)rn  =  /  E(x,  t2) 

datX  ^  I'  )  D,m  ^ 

^(h)m  =  f  F{x,  h), . /(*)„, m*i 

a(ti)„  =  ACj((i)a„  +  SC„j(ti)#ll 

/?(*i)n  =  ADn9(h)D,n  +  BDng(ti)Ntn 

7fe)n  =  ACng{t2)Dn  +  BCng(t2)Nn 
A{t2)n  =  ADng(t2)Dn  +  BDng(t2)Nn 


f  an  Jn\  =  (  ACn  BCn  \  g(t i)Dn  g(h)Dn  \ 
\  fin  ^n  )  V  ADn  BDn  J  9^1)  N,n  g{t 2)jv,n  / 

f  ACn  BCn  \  =  (  Oin  \  (  g(t i)A„  9(h)  DiU  ^ 

\  ADn  BDn  J  fin  An  J  y  9(h)  NjU  g(h) jv,n  / 
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Figure  1:  Guitar  String 


Figure  1:  Dual  Time  Problem  -  similar  to  the  Initial  Value  Problem 
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Dual  Position  Problem  (conjugate  to  time  problem) 


F(x,t)o,n  =  f{x)ng(t)n 

oo  r  i  r  i 

F(x -  =  S  A„g{t)  +  Bng(t)  CJ{x)  +  DJ{x) 

n= 0 


OO  r 

F(X’t)data  —  ACngDnf Dn BCngNnf Dn  +  ADngDnf Nn  +  BDngNnf Nn 


a(xi)m  =  I  F(Xl,t)datag(t)Dmdt 


„  oo  r  i 

a(Xl)m  —  J  ACngDJD,n  +  BCngNJD,n  +  ADn9DJN,n  +  BDn9NJN,n  9( 

n  0 

a(xi)„  =  ACJ(xi)Dn  +  ADnf(xi)Nn 


a(Xl)m  =  j  F{Xl,t)iatag{t)Dmdt 
P(.xi)m  =  f  F(Xl,t)datag(t)Nmdt 

7  {x2)m  =  j  F(x2,t)datag(t)Dmdt 

A(x2)m  =  J  F(X2,t)Mj(t)Nmdt 
<x(xi)n  =  ACng(t\)Dn  +  ADng(ti)N7 

P{xi)n  =  BCng{t1)Dn  +  BDng(tx)N 

l{x2)n  =  ACng(t2)Dn  +  ADng(t2)Nt 

A(x2)n  =  BCng{t2)Dn  +  BDng(t2)N 


(an  7„\  =  [ACn  ADn\  f(Xl)Dn  f(xt)D\ 
vA  An)  \BCn  BDnj  f{x-i)Nn  f(x2)Kn) 

I  ACn  ADn\  =  7„\  /(*i)D,„  /(®2)A„V 

v  BC'n  BDn  J  (  &  AnJ  (  /(^i)jv,„  /(®2)jv,„  j 
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Figure  3:  Dual  Position  Problem  -  conjugate  to  the  Dual  Time  Problem 


Figure  4:  Having  found  the  amplitudes,  the  solution  is  projected  across  the  spatial  domain. 
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Multiple  Position  Problem 


F(x,t)0,n  =  f{x)ng(t)n 

OO 

F(X^)data  =  ^  N,n  ~ ^ nt  ~ "^n]  \C nf  j^  n  n  <f  N  ,n  C*  r 


u\j 

data  AC nQ  Dnf D,n  ^  C nQ  Nn  f D  ,n  ~\~  AD  ^  D n/ ^ n  +  B  D  ^  N  n  f N  n 


F{x,t)data  =  •••  16  terms ,  needs  8  locations 
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X 


X 

8 


5:  Having  found  the  amplitudes,  the  solution  is  projected  across  the  spatial  domain. 
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a  for  mode  m=1 5  out  of  50  modes 
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Dirichlet  Mode  35  (350X1185) 
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Neumann  Mode  4  (350X1185) 
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Neumann  Mode  100  (350X1185) 
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Methods  of  solution 

Zel’dovich  (1985):  Velocity  vectors  fields  can  be 
extracted  from  two  scalar  potentials 


u 

=  Vx 

1 

L  v 

V2 

n  ^  n 

V2 

m  1  m 

x¥d  I  =  0 

woundary 

(« •  V® " )  =  0 


Lipphardt  et  al.  (2000):  Addition  of  forcing  terms 
allows  for  non-conservation  of  mass  through  a 
boundary  -  ie.  Water  from  rivers  or  the  ocean  is 
accounted. 

V  20  (x,  y,  0,  t)  =  S  @  (t) 

(jfl  VO)  |  boundary  (jTl  ^m0(j  t>/ )  I  boundary 


Putting  It  All  Together 

•  'Fis  the  stream  potential  (vorticity  mode). 

•  O  is  the  velocity  potential  (divergent  mode). 

•  Situation  analogous  to  (E,B)  fields  from  E&M. 

•  The  vector  field  representation  can  be  separated 
into  two  eigenvalue  equations. 

•  Source  term  solved  via  Poisson’s  equation. 

•  The  total  vector  field  is  written  as  a  sum  over  all 
states  for  each  representation. 

N  M 

(W,V)  =  2X  (Un,Vn)D  +  Yjbm(Um’VJN  +  («(0 Mt))s 

n=l  m= 1 


Time  Series  Analysis 

•  Chesapeake  flow  can  be  written  as  a 
Normal  Mode  expansion  =>  u(r,t) 

•  N 

Mo  =  +  bn(t)untN(r)+  u(r,t)src 

n=  1 

•  Use  Galerkin  method  to  extract  a(t),  b(t) 

•  Due  to  limitations  in  data  collection: 

-  use  QUODDY  (model  based  on  data) 

-  Partial  domain  Galerkin  method 


Partial  Domain  Galerkin  Method 

f(r,t)=  5X(0'*'„  Dif)  +  bn(t)®nN(r) 

n- 1 

^  C'1)  •  ^  (^>  + 

S„  =2>„  (0  { J  (?)  ■ ■ 

72=1  72=1 

let 

/L,  =  j(l \{r)-^m{r)  d£lpartiai 

N  N  ■ 

am  =  +  Z  {^m}^(0 

72=1  72=1 


Histogram  of  Orthonormalization  of  T.  vs.  T. 
at  Highest  Resolution 


1  /  N  J  T.  x  'P.  c/x  dy 


Fig.  6  Orthogonality  between  any  two  modes  shows  the  degree 
which  information  about  the  system  overlaps.  For  a  basis  set,  the 
requirement  is  that  all  functions  spanning  the  space  have  a  zero  in¬ 
tegral  when  considering  the  product  of  any  two  modes  over  the  do¬ 
main.  The  Neumann  modes  for  the  Bay  are  decent  as  compared  to 
similar  calculations  on  the  unit  square  and  circle. 


Employing  the  same  techniques  for  the  Neumann  modes, 
agreement  could  not  be  reached  within  an  acceptable  range 
(10-20%).  While  studying  the  unit  square  and  circle  in 
order  to  establish  baseline  performance,  the  finite  differ¬ 
ence  results  typically  varied  more  than  FEMLAB,  so  most 
likely  the  error  is  coming  from  the  finite  difference  scheme 
and  not  FEMLAB.  Although  not  conclusive,  the  agree¬ 
ment  between  the  two  methods  indicates  stability  of  the 
solution  set. 

Although  not  shown  in  this  report,  the  Dirichlet  and 
source  terms  were  calculated  for  the  Chesapeake  Bay  and 
can  be  found  in  references  [8]  and  [11].  Comparison  with 
the  finite  difference  scheme  will  require  further  study  to 
validate  either  the  FEMLAB  result  or  the  finite  differences 
result.  Given  that  Neumann  conditions  are  better  suited  for 
finite  element  analysis  leads  one  to  trust  the  FEMLAB  re¬ 
sults,  however,  improvements  to  the  finite  differences  are 
still  possible,  which  will  help  validate  the  FEMLAB  so¬ 
lution  set.  The  solutions  arising  from  the  source  term  be¬ 
haved  as  expected.  For  a  detailed  inquiry  into  the  complete 
eigenmode  set  calculated  to  date  for  the  Chesapeake  Bay, 
consult  reference  [8].  These  results  can  also  be  seen  at  the 
website  listed  at  USNA  [9].  Future  work  on  the  velocity 
vector  field  for  the  Chesapeake  Bay  will  include  full  three 
dimensional  analysis  using  FEMLAB  and  Quoddy  in  con¬ 
junction. 

Concurrent  with  this  study,  the  residence  time  of  par- 
ticulants  was  analyzed  using  the  Navier-Stokes  equations 


Eigenvalues  of  Chesapeake  Bay  Neumann  Modes 


Fig.  7  The  progression  of  eigenvalues  for  the  Neumann  modes  of 
the  Chesapeake  Bay.  No  fluctuation  is  observed  as  the  mode  number 
reaches  100,  indicating  that  the  feature  size  of  these  modes  is  far 
from  the  size  of  the  mesh. 

and  integrating  the  velocity  fields  [10].  By  comparing  trends 
in  velocity  fields  from  an  eigenmode  calculation  to  a  Navier- 
Stokes  calculation,  one  can  analyze  those  parts  of  the  ve¬ 
locity  domain  that  give  rise  to  certain  behaviors.  This  com¬ 
parison  represents  the  initial  steps  needed  to  perform  a 
full  Normal  Mode  Analysis.  Plans  have  begun  to  further 
instrument  the  Bay,  leading  to  the  question,  how  many 
real-time  measurements  are  required  in  order  to  predict 
important  behavior  within  the  Chesapeake.  Normal  Mode 
Analysis  will  be  used  to  help  answer  this  question. 
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Fig.  18  The  mesh  from  QUODDY  [9].  The  highest  density  mesh  is  Fig.  19  The  mesh  from  QUODDY’s  boundary  remeshed  within 

in  the  region  of  greatest  change  in  the  depth.  COMSOL.  The  highest  density  mesh  is  in  the  region  of  greatest 

change  in  curvature  near  the  boundary. 


cussion  of  fractals!  The  boundary  for  the  calculations  for 
this  paper  were  taken  from  the  QUODDY  model  that  the 
meshes  will  match.  Clearly,  the  feature  set  of  the  boundary 
has  been  reduced  by  averaging  the  locations  of  the  edges, 
effectively  smoothing  the  Chesapeake.  Even  so,  the  min¬ 
imal  mesh  for  COMSOL  is  3200  elements.  By  adding  in 
features  from  the  QUODDY  mesh,  the  coarsest  mesh  will 
easily  be  over  5000  elements. 

A  quick  review  of  oceanographic  websites  indicates 
that  there  are  10  data  monitoring  stations  taking  current 
data  in  the  Chesapeake.  This  sparse  amount  of  data  sug¬ 
gests  a  data  coverage  of  approximately  0.2%.  One  pos¬ 
sibility  for  increasing  the  data  coverage  is  to  have  many 
more  data  taking  stations.  This  is  costly  and  impractical 
as  the  Chesapeake  is  a  major  waterway  for  commerce  and 
military  use.  An  alternative  approach  is  to  take  a  truly  poor 
mesh  of  the  Chesapeake  Bay.  By  severely  lowering  the 
mesh  resolution,  the  data  coverage  will  increase  at  the  ex¬ 
pense  of  numerical  accuracy.  By  time  averaging  data,  the 


longer  time  taken  can  also  overcome  poor  spatial  cover¬ 
age.  Each  of  these  effects  suggest  an  uncertainty  princi¬ 
ple  8 (numerics)  8 (spatial) 8 (temporal)  —  constant.  The 
standard  checks  for  such  an  analysis  will  be  limited  to  sim¬ 
ple  dimensional  analysis  or  1-D  signal  applications,  as  in 
section  4. 

COMSOL  is  taking  projects  such  as  the  Chesapeake 
Bay  into  new  territory  as  the  complexity  rises.  The  bank 
of  tests  from  simple  systems  may  not  easily  apply  to  the 
richness  of  the  systems  capable  of  being  calculated  by 
COMSOL,  yet  difficult  to  analyze  when  combined  with 
real-world  data.  This  is  an  on-going  project.  Stay  tuned. 
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Conclusion 

•  Possibly  found  way  to  use  ~10  monitoring  stations  to  extract 
full  Chesapeake  Bay  flow  field. 

•  Time  series  has  a  good  chance  to  work! 

•  Due  to  COMSOL  and  usage  of  fuller  geometric  solutions, 
challenge  older  precepts  based  on  signal  processing. 

•  If  orthonormality  is  the  strongest  requirement,  create  post¬ 
processing  options  to  massage  results  obtained  from 
solvers.  Given  them  tolerances  to  adjust  results. 

•  Please  allow  for  easier  adjustment  and  creation  of 
meshes! 

•  Visit  us  at:  http://web.usna.navy.mil/~rmm/ 


Thanks 
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